A fast algorithm for computing minimal-norm solutions to underdetermined systems of linear equations

We introduce a randomized algorithm for computing the minimal-norm solution to an underdetermined system of linear equations. Given an arbitrary full-rank matrix Am×n with m < n, any vector bm×1, and any positive real number ε less than 1, the procedure computes a vector xn×1 approximating to relative precision ε or better the vector pn×1 of minimal Euclidean norm satisfying Am×n pn×1 = bm×1. The algorithm typically requires O(mn log( √ n/ε)+m3) floating-point operations, generally less than the O(m2 n) required by the classical schemes based on QR-decompositions or bidiagonalization. We present several numerical examples illustrating the performance of the algorithm.